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ABSTRACT 

We study the dynamics of phase transitions in the interstellar medium by means of three-dimensional hydrody¬ 
namic numerical simulations. We use a realistic cooling function and generic nonequilibrium initial conditions to 
follow the formation history of a multiphase medium in detail in the absence of gravity. We outline a number of 
qualitatively distinct stages of this process, including a linear isobaric evolution, transition to an isochoric regime, 
formation of filaments and voids (also known as “thermal” pancakes), the development and decay of supersonic 
turbulence, an approach to pressure equilibrium, and final relaxation of the multiphase medium. We find that 
l%-2% of the initial thermal energy is converted into gas motions in one cooling time. The velocity field then 
randomizes into turbulence that decays on a dynamical timescale oc t~ a , 1 < a <2. While not all initial con¬ 
ditions yield a stable two-phase medium, we examine such a case in detail. We find that the two phases are well 
mixed with the cold clouds possessing a fine-grained structure near our numerical resolution limit. The amount 
of gas in the intermediate unstable phase roughly tracks the rms turbulent Mach number, peaking at 25% when 
M. nm ~ 8, decreasing to 11% when M. rms ~ 0.4. 

Subject headings: hydrodynamics — instabilities — ISM: structure — turbulence 


1. INTRODUCTION 

Thermal instability (TI) has many implications in astro¬ 
physics (e.g., a clumpy interstellar medium [ISM], stellar atmo¬ 
spheres, star formation, globular cluster and galaxy formation, 
etc.; see Meerson 1996 for a recent review). The instability may 
be driven by radiative cooling of optically thin gas (radiation- 
driven TI) or by exothermic nuclear reactions (Schwarzschild 
& Harm 1965). 

Linear stability theory for a medium with volumetric sources 
and sinks of energy in thermal equilibrium was developed by 
Field (1965), who identified three unstable modes: the iso¬ 
baric mode (the pressure-driven formation of condensations not 
involving gravitation) and the two isentropic modes (the over¬ 
stability of acoustic waves propagating in opposite directions). 
Hunter (1970, 1971) extended these results to an arbitrary non¬ 
stationary background flow, showing that cooling-dominated 
media are potentially more unstable than those in equilibrium, 
while heating provides stabilization. The most common appli¬ 
cations of thermal instability to the ISM and star formation deal 
with the isobaric mode that was employed to explain the ob¬ 
served multiphase structure of the ISM (Pikel’ner 1968; Field, 
Goldsmith, & Habing 1969; McKee & Ostriker 1977; McKee 
1990; Wolfire et al. 1995). 

Analysis of infinitesimal perturbations gives two characteris¬ 
tic length scales for the isobaric condensation mode: (1) a cool¬ 
ing scale A p = c / uj v (where c is the adiabatic sound speed and w p 
is the growth rate) and (2) a critical scale X K = c\Jtd/oJ v (where 
td is the characteristic thermal diffusion time). These two length 
scales define short-, intermediate-, and long-wavelength lim¬ 
its (Meerson 1996; Kovalenko & Shchekinov 1999). In the 
short-wavelength limit, small isobaric perturbations are inhib¬ 
ited by heat conduction, so that w p < 0 for A < X K . In the 
long-wavelength limit, the perturbations cannot grow isobari- 
cally because of the finite sound speed effects, and thus u; p —> 0 
for A/A p —* oo. This is only true if the gas is isochorically sta¬ 
ble (Parker 1953; Field 1965; Shchekinov 1978), otherwise 


the growth rate remains finite: w p —► ui p > 0 for A/A p —> oo, but 
only large-scale temperature perturbations are growing, thus re¬ 
sulting in pressure variations and the formation of shock waves. 
The growth rates and characteristic scales depend on the heating 
and cooling properties of a given medium. Under the ISM con¬ 
ditions, if one assumes thermal equilibrium (i.e., an exact bal¬ 
ance between cooling and heating), isochoric instability man¬ 
ifests itself only at relatively high temperatures, T > 10 5 K. 
However, in cooling-dominated regimes, it can develop at tem¬ 
peratures as low as 10 3 K. 

A specific feature of TI in the ISM is a large ([> 1 dex) gap in 
gas densities between the two stable phases. The density range 
of interest in a galaxy formation context is even larger. This 
implies the importance of nonlinear effects in the dynamics of 
phase transitions. Nonlinearity brings into play nonequilibrium 
effects. Already weakly nonlinear development of condensa¬ 
tions in an initially homogeneous gas in thermal equilibrium 
drives the system away from equilibrium. The mean pressure 
drops since p 2 > p 2 and cooling overcomes heating globally 
(Kritsuk 1985). Later, on a timescale of ~ w p , as condensa¬ 
tions get denser and cooler, the isobaric condition A <C A p be¬ 
comes violated locally within them, so the system departs from 
pressure equilibrium. These effects are essential for the isobaric 
mode of TI in the ISM and star formation contexts. There¬ 
fore, analytical nonlinear solutions to “isobarically” reduced TI 
equations are insufficient to describe the radiative stage of the 
phase transition (f ~ w” 1 ). 

During this strongly nonlinear stage large-scale condensa¬ 
tions form in such a way that gas moves almost inertially and 
its kinetic energy dominates thermal energy (p <C pv 2 ). Ac¬ 
cordingly, the gas velocities in these condensations are of the 
order of the sound speed in the unperturbed state. The situ¬ 
ation here is directly analogous to the long-wave gravitational 
instability, so that results concerning the origin of cellular struc¬ 
ture and Zel’dovich (1970) “pancakes” can be entirely carried 
over to the case of long-wave TI (Meerson & Sasorov 1987). 
Sasorov (1988) gave an elegant proof that qualitatively the 
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FlO. 1.— Snapshots of the gas density field (perspective volume rendering): (a) First condensation at t — 0.07 Myr, (b) thermal pancakes at / = 0.1 Myr, (c) 
collapse and turbulization of cellular structure at t = 0.17 Myr, (d) two-phase medium at f = 1.5 Myr (5 pc box, 256 3 grid points). The log density color coding is as 
follows: The dense blobs at the intersections of the filaments, p > 10 -2- g cm -3 , are light blue; the stable cold phase, p £ [ 10 23 ,10 32 ] g cm -3 , is blue; the unstable 
density regime, p £ [ 10 23 7 ,10 3311 1 g cm -3 , is yellow to brown; and the low-density gas, including stable warm phase (p < 10 33 7 g cm -3 ), is a transparent red. 
The figure is also available as an mpeg animation in the electronic edition of the Astrophysical Journal. 


same result applies to small-scale TIs; i.e., the onset of TI pro¬ 
duces voids and highly flattened condensations along certain 
two-dimensional surfaces. These are also called “thermal” pan¬ 
cakes (Meerson 1996). The formation of filaments was simul¬ 
taneously noticed in two-dimensional numerical simulations of 
TI in the solar transition region (Dahlburg et al. 1987; Karpen, 
Picone & Dahlburg 1988). Ever since, thermal pancakes are be¬ 
ing rediscovered both analytically and in numerical simulations 
(e.g., Lynden-Bell & Tout 2001). 

Thermal pancakes are transient. However, what happens 
next, before the evolution turns to a conductive relaxation stage 
(Meerson 1996), until recently has remained the “terra incog¬ 
nito” of TI theory. The problem of “postradiative” mechanical 
relaxation toward a static multiphase medium requires a solu¬ 
tion for the full set of hydrodynamic equations that can only be 
obtained numerically. One-dimensional simulations pioneered 
by Goldsmith (1970) demonstrated that TI develops large mo¬ 
tions in the ISM (see also Hennebelle & Perault 1999 for an 
example of how large motions can trigger TI). For some time, 
progress in this direction was precluded by numerical diffi¬ 
culties in modeling convergent cooling flows with high Mach 
numbers and high-density contrasts (e.g., Vazquez-Semadeni, 
Gazol, & Scalo 2000). Recent multidimensional numerical 
simulations of the ISM evolution in disklike galaxies include 
effects of gravity, differential rotation, star formation and su¬ 
pernova feedback (de Avillez 2000; Vazquez-Semadeni et al. 
2000; Wada & Norman 2001; Wada 2001). However, it is hard 
to determine the role of TI in shaping the ISM structures found 
in these models, partly because the simulations still do not re¬ 
solve length scales important for TI and partly because of the 
additional physical effects. 

The purpose of this Letter is to report on results of three- 
dimensional numerical simulations of classical TI that fill the 
gap in theory, exploring in detail the late radiative stage and 
postradiative relaxation toward a multiphase medium. Our ma¬ 
jor result is that formation of thermal pancakes induces turbu¬ 
lence in the ISM that serves as a nonlinear saturation mecha¬ 
nism for TI. As a consequence of a turbulent cascade, (1) infor¬ 
mation about initial perturbations is lost, including the imprints 
of heat conduction in the density power spectrum during the 
linear stage, and (2) turbulent diffusion becomes the dominant 
transport mechanism during the postradiative relaxation stage. 


2. SIMULATIONS 

We solve the equations of ideal gasdynamics (eqs. [6]- 
[9] in Field 1965) in a cubic domain with periodic boundary 
conditions, assuming zero conductivity and no gravity. The 
generalized cooling function £ contains cooling and heating 
terms, £ = pA(T,Z) — T. The temperature dependence of ra¬ 
diative cooling A (T,Z) for T £ [10,10 8 ] K is provided by M. 
Spaans (2000, private communication; see also Wada & Nor¬ 
man 2001). The gas metallicity Z is assumed to be solar in 
most of the models, but we also considered a case with a subso¬ 
lar value of O.IZq. For simplicity we assume that the heat¬ 
ing rate I' = const, and we define it in units of the cooling 
rate in an initial unperturbed state: T = QpqA{Tq,Z), where 
Q £ {0,0.01,0.1,0.5,1}. The initial unperturbed state is uni¬ 
form with gas density po £ {0.167,1.67} x 10“ 24 g cm -3 and 
gas temperature Tq £ {7 000, 2 x 10 6 } K. The box size L £ 
{5,100,500} pc was chosen such that at least the cooling scale 
A p ,o for initial state is well resolved on a discrete grid of 128 3 or 
256 3 cells. In most cases, we initially impose isobaric density 
perturbations with a three-dimensional power spectrum index 
of -3, an rms deviation of e = 0.05, and a high-wavenumber 
cutoff k max £ {8,32},' assuming zero initial velocities. A few 
cases were started up with random velocity perturbations, vary¬ 
ing the amplitude and assuming uniform initial density. 

The simulations were performed using two hydrosolvers im¬ 
plemented in the Enzo code developed by Bryan & Norman 
(2000): one is a direct Eulerian version of the PPM scheme 
(Colella & Woodward 1984), and the other is the finite differ¬ 
ence Eulerian method originally implemented in ZEUS (Stone 
& Norman 1992). 


3. RESULTS 

Here we present results for our fiducial model computed on 
a “moderate”-resolution grid of 256 3 cells assuming /, = 5 pc, 
po= 1.67 x 10 -24 gcm -3 , T ( ) = 2x 10 6 K, Q = 0.1, solar metallic¬ 
ity, and isobaric perturbations with e = 0.05, k max = 32 (Fig. 1). 
With these initial conditions, the linear isobaric mode develops 
on a timescale au} « 0.3 Myr, the cooling length A Pi o ~ 70 pc, 
and the critical length for heat conduction A Ki o ~ 8 pc. Thus, 
our initial perturbations formally drop into a short-wavelength 


1 These are normalized values; k max = 32 for a grid 256 3 implies that the cutoff scale spans eight grid zones. This is slightly below the dissipation limit of our 
numerical scheme. 
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Fig. 2.— Snapshots of phase diagrams (timing and labels correspond to those in Fig. 1): (a) t = 0.07 Myr, (c) t = 0.17 Myr, (e) t = 0.5 Myr, (d) t = 1.5 Myr. The 
black dots show scatter plots of pressure vs. density. The “dash” at P = 4.55 x 10 -10 dyn cm -2 in (a) shows the isobaric initial conditions. Background yellow-filled 
contours specify the part of the phase plane where isobaric mode is unstable; overlaid magenta contours are the regions of isochoric instability. The thick solid line 
shows thermal equilibrium curve. Density PDFs are plotted at the bottom of each panel (see scale to the right). 


limit . 2 Assumed cooling and heating rates and the value of 
the initial gas density imply the existence of attracting bistable 
thermal equilibrium, with physical conditions close to those in 
the ISM (Fig. 2), although the pressure range for bistability is 
somewhat narrower than in the “standard” model of Wolfire et 
al. (1995). 

With our setup, heating does not compensate cooling in the 
initial state; thus, TI sets in on a time-dependent background 
state. As the gas cools, the instability channels a part of its ther¬ 
mal energy into kinetic energy of converging flows that create 
condensations and evacuate gas from underdense voids (Fig.3). 
Growing density variance accelerates cooling efficiency, so the 
gas as a whole cools much faster than a homogeneous medium 
with the same initial setup. Since initial perturbations span a 
range of linear scales, the epoch of the thermal pancakes’ for¬ 
mation lasts from ~ 0.07 Myr, when the first cold blob forms, 
to ~ 0.1 Myr, when the mean kinetic energy and then pres¬ 
sure variance reach their maximal values (see Figs, la, lb, 2a, 
and 3). 

The pancakes themselves exhibit a rather complex inner 
structure. As the isobaric compression gives way to isochoric 
cooling (Kritsuk 1990; Burkert & Lin 2000), accretion shocks 
develop within the cooling condensations at temperatures cor¬ 
responding to two strips of stability in the phase plane (Fig. 
2a). Dense gas in the cores of condensations cools further as 
it contracts in a regime similar to the explosive condensation 
described in Meerson & Sasorov (1987); see Figure 2a. 

Due to asymmetries in the initial conditions, the dense cores 
gain nonzero momentum that drives a bottom-up collapse of the 
hierarchical cellular structure composed of thermal pancakes. 
As a result, a single large void forms in the periodic box, and 
p m i n attains its global minimum by / = 0.13 Myr (Figs, lc and 
3). At this point, the gas density variance approaches its maxi¬ 


mum, and the density spans a range of 5.5 dex. The highest den¬ 
sity gas quickly relaxes to thermal equilibrium. However, the 
dense cold blobs reexpand slightly until they reestablish pres¬ 
sure balance with the less dense environment, so the mass frac¬ 
tion of the gas in cold stable phase “H” decreases 3 (see Fig. 3). 

By t ~ 0.15 Myr, the information about the details of the 
initial perturbations is lost, and highly compressible supersonic 
turbulence with an rms Mach number of about 10 is fully de¬ 
veloped in the computational domain (Fig. 3). The velocity 
power spectrum fills in at small scales as the turbulent cascade 
settles and then only slightly evolves on a time scale of ~ 1 Myr 
as turbulence decays [£* ss 9.0 x 10 13 (//0.3 Myr) 2 ergs cm 3 
fort £ (0.3,1.7) Myr, cf. Mac Low, Klessen, & Burkert 1998]. 
The density probability distribution function (PDF) exhibits a 
power-law excess at high densities, where the “effective” equa¬ 
tion of state is soft 4 . Low-density gas is nearly adiabatic (7 = |); 
therefore, the PDF also has a power-law regime at low p (Figs. 
2c and 2e). The power-law regimes appear in response to devi¬ 
ations from a special isothermal case in which lognormal PDFs 
are produced due to the fact that the local Mach number is in¬ 
dependent of the density (Passot & Vazquez-Semadeni 1998). 

Relaxation timescales toward thermal and dynamical equilib¬ 
ria depend on local physical conditions. Cold dense gas quickly 
settles to a thermal equilibrium, but its dynamical relaxation 
proceeds quite slowly (by t = 1.5 Myr, its kinetic energy domi¬ 
nates the thermal energy by up to a factor of 3 in the most dense 
blobs). Warm tenuous gas, instead, undergoes a fast transition 
to a subsonic regime with a quite uniform density distribution 
(see PDFs in Figs. 2e and 2d), but it takes longer to establish 
thermal equilibrium. As soon as the system relaxes to a bistable 


2 Since initial perturbations are not infinitesimal, —\max,a = 0.23. TI reaches the strongly nonlinear regime by 0.07 Myr. For such perturbations, both A Pi o and X K p 
are effectively lower than their linear estimates. 

3 We follow the notation of Field et al. (1969). 

l_j_ 

4 For A (T) oc T a and T = const, the equilibrium pressure p eq oc p eq a , and since a > 1 at low temperatures, the effective adiabatic index j eq < 1. 
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Fig. 3.— Time evolution of global variables. Top left panel: p max (blue), p m i n (red), < p 2 > / < p > 2 (black), 10 < p 2 > / < p > 2 (magenta, scale is at the right). 
Middle left panel: Total energy (black), thermal energy (red), kinetic energy (blue). Bottom left panel: Mass-weighted rms Mach number (blue), rms enstrophy 
(green). Top right panel: Mass fractions of different thermal phases: hot (labeled HIM, T > 19,000 K) warm stable (labeled F, T £ [8000,19,000] K), intermediate 
unstable (labeled G, T E [600,8 000] K), and cold stable (labeled H, T < 600 K). Bottom right panel: Same as above, but for volume fractions. 


state with pressure close to the allowed minimum, 5 the density 
PDF shows a bimodal distribution, phase fractions by mass sat¬ 
urate at ~ 42% for the warm phase, and ~ 44% for the cold 
phase, and ~ 14% of the mass falls in the unstable temperature 
regime. The substantial amount of gas in the unstable regime 
can be explained as being due to turbulent diffusion. Vorticity 
is generated by baroclinic instabilities and in shocks, which de¬ 
velop during the radiative stage at the interface of growing con¬ 
densations (Dahlburg et al. 1987; Karpen et al. 1988). The val¬ 
ues we found for the unstable gas fraction are somewhat lower 
than those obtained in two-dimensional models by Gazol et al. 
(2001). To see whether turbulence was at least partially driven 
by mass exchange between phases, we zeroed the velocity field 
at t ~ 1.5 Myr. Residual pressure variations regenerated the tur¬ 
bulence, but at a much lower level. The gas mass fraction in the 
unstable regime “G” shrunk from 14% to only 2%; the remain¬ 
der is incorporated into the cold stable phase H. We conclude 
that Tl-induced turbulence is purely decaying in our simula¬ 
tions. 

Since at least some part of the gas does not stay in thermal 
equilibrium (Fig. 2d), the exact fraction of thermally unsta¬ 
ble gas should be determined by Hunter’s (1970) TI criterion. 
A rough estimate, however, could be obtained using density 
and temperature limits for the stability of thermal equilibrium, 
which are used in Figs 1 and 3, respectively. Note that such es¬ 
timates are robust for only a medium in a quasi-isobaric state. 
One has to be cautious obtaining temperature estimates for the 
unstable gas based on the assumption of thermal equilibrium 
(cf. Heiles 2001a). 

Only because we were able to resolve the length scales that 
were smaller than the cooling scale A p in the vicinity of the 
bistable regime did the gas in the box relax to an isobaric 
state. 6 Seed perturbations in the corresponding range of lin¬ 
ear scales are always available to feed the instability, provided 
the turbulent cascade exists. The resulting two-phase medium 
is dynamic. Turbulent velocities do not correlate with density. 


and dense clouds appear to be shapeless random aggregations 
of cold Lagrangian gas parcels, forming a “fractal” substrate 
(Fig. Id). ' 

4. DISCUSSION 

Our fiducial case was constructed to produce a stable two- 
phase medium because of its relevance to the Galactic ISM. We 
are interested in TI over a wide range of conditions as might be 
found in the ISM of high-redshift protogalaxies. We have sim¬ 
ulated other cases with different parameter choices that do not 
produce stable two-phase media. However, we find they all de¬ 
velop turbulence in the nonlinear radiative stage of TI. Here we 
briefly discuss how the turbulence and asymptotic phase struc¬ 
ture depend on initial conditions, deferring a more complete 
discussion to a future paper (A. Kritsuk & M. Norman 2002, in 
preparation). 

The level of induced turbulence is determined by the effi¬ 
ciency of conversion of the initial thermal gas energy into ki¬ 
netic energy of turbulent flow by nonlinear development of TI: 
E" ax = C(po,T( h e.Q,L)E th (0) (see Fig. 3). The conversion fac¬ 
tor C is a complex function of its variables. It varies from about 
2% to 1% in our models. In general, higher e and/or To val¬ 
ues provide higher conversion; a lower heating level (Q < 1) 
supports TI and therefore works in the same direction. Larger 
boxes, as a rule, also produce more turbulence. The turbulence 
is induced on the initial cooling time and decays on a dynam¬ 
ical timescale, which is typically much longer. The turbulent 
Mach number at t ~ w p depends on the mean temperature 
at this time. For nonequilibrium initial conditions {Q <C 1 or 
Q = 0), this is much lower than the initial temperature. In our 
fiducial case, M rms peaks at 8, dropping to 0.4 after 20 initial 
cooling times. For equilibrium initial conditions (Q = 1), tur¬ 
bulent velocities remain subsonic (Ai nm ~ 0.3). But we would 
expect higher Mach numbers if the bistable range of pressure 
were wider than provided by our adopted cooling function. 

The initial gas density determines the number and mass frac- 


5 The pressure hovers near P,„,„ because of the shape of the thermal equilibrium curve, which is in turn controlled by our assumed heating and cooling functions. 

6 It was the necessity to resolve isobaric mode in the vicinity of bistable thermal equilibrium that forced us to choose the box size L = 5 pc. 
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tions of thermal phases in the relaxed state depending on its 
position relative to the valleys and hills on the thermal equi¬ 
librium curve. This is consequence of our choice of constant 
volume boundary conditions, which means that the mean den¬ 
sity in the box remains constant. After the rapid cooling stage, 
our models with low initial densities po = 1-5 x 10 -25 g cm -3 , 
high temperature Tq = 2 x 10 6 K, and Q G {0.3,1} generate 
turbulence, evolve through a transient three-phase stage, and 
then relax to a single-phase low-pressure warm ISM. While 
turbulence is a generic feature of nonlinear saturation of TI, 
our simulations show that detailed turbulent properties and the 
nature of emerging multiphase medium do depend sensitively 
on the Mach number and effective equation of state controlled 
by heating and cooling; this will be discussed elsewhere. Two 
identical simulations, except that cutoffs in initial power spec¬ 
tra were different (k max = 8 and 32 on a 128 3 grid, L = 100 pc, 
Q = 0), demonstrated considerable structural differences in den¬ 
sity distributions at the thermal pancake stage, t lp , and surpris¬ 
ingly similar “chaotic” density structures and identical veloc¬ 
ity power spectra at ~ 6t tp , when turbulent mixing covered the 
whole computational domain. This implies that the imprints of 
heat conduction in the density power spectrum during the linear 
stage could be erased later by the developing turbulent cascade. 

TI is certainly not the only potential source of turbulence in 
the ISM, but it cannot be ignored at least in those scenarios that 


actively employ TI to explain the origin and properties of ob¬ 
served objects. We suggest a paradigm shift concerning the role 
of thermal instability in the ISM and the nature of multiphase 
ISM. The idea of “static” two-phase ISM introduced in late 
1960s (pressure-confined thermally stable dense clouds embed¬ 
ded in rarefied intercloud gas forming as a result of TI and sub¬ 
ject to phase exchange due to cloud evaporation/condensation) 
must give way to the notion of a dynamic multiphase ISM, in 
which TI induces slowly decaying turbulence and in which tur¬ 
bulent diffusion regulates phase exchange processes. In this 
new emerging picture, the dense clouds are shapeless random 
aggregations of cold Lagrangian gas parcels; the clouds do 
not preserve their identity in real space on their sound cross¬ 
ing timescale until self-gravity tightens the fragments up into a 
self-gravitating cloud to form stars. Our results may suggests 
modifications to the scenario of a three-phase ISM (McKee & 
Ostriker 1977; McKee 1990; Heiles 2001b) that are yet to be 
understood. 

We are grateful to Marco Spaans for providing cooling func¬ 
tions prior to publication. We acknowledge useful conversa¬ 
tions with George Field and Chris McKee. This work was par¬ 
tially supported by NRAC computer grant MCA98N020N and 
utilized the NCSA Silicon Graphics Origin2000 system at the 
University of Illinois at Urbana-Champaign. 
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